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1.  Introduction 
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Experiments  with  airborne  gravimeters  have  been  performed 
over  the  last  15  years.  Although  the  equipment  designed  for 
this  objective  has  a high  degree  of  sophi sti cati on, the  results 
obtained  so  far  are  not  accurate  enough  for  geodetic  purposes. 

The  reason  lies  in  the  complicated  structure  of  the  force  field 
acting  on  a moving  gravimeter.  Gravimeters  are  basically 
accelerometers  and  they  measure  the  resultant  of  gravitational 
and  inertial  forces.  If  they  are  used  as  stationary  instru- 
ments, as  in  most  terrestrial  applications,  the  only  inertial 
force  acting  on  the  gravimeter  is  the  centrifugal  force.  There- 
fore, the  output  of  the  instrument  is  the  combined  effect  of 
gravitational  attraction  and  centrifugal  force,  i.e.  gravity. 

The  situation  is  more  complicated  in  a moving  gravimeter.  The 
Coriolis  force  has  to  be  taken  into  account  and,  more  important, 
irregular  accelerations  of  the  base  will  strongly  influence  the 
result  of  the  measurements.  Such  undesirable  inertial  forces 
are  especially  strong  in  a moving  aircraft  and  there  is  no  way 
to  rigorously  separate  the  gravitational  part  from  the  inertial 
part  by  using  gravimeter  measurements  only.  Therefore,  additional 
information  is  necessary  to  extract  the  gravitational  effect. 

Two  methods  have  been  proposed  to  reach  this  goal.  In 
the  first  one  information  on  the  frequency  behaviour  of  the 
different  forces  is  used  to  separate  gravity  and  disturbing 
accelerations  by  statistical  filtering  techniques.  Meissl 
(1970)  has  investigated  this  approach  using  the  theory  of  sto- 
chastic processes  and  certain  assumptions  on  the  power  spectra 
of  the  force  fields  involved.  He  concludes  that  it  is  most  diffi- 
cult to  separate  gravitation  and  inertia  in  the  medium  fre- 
quency range  with  half  wavelength  between  30  and  150  km.  This 
is  only  possible  if  detailed  information  on  the  two  spectra  is 
available  which  usually  will  not  be  the  case.  The  high  fre- 
quencies can  be  blocked  by  a low  pass  filter,  the  low  fre- 
quencies can  be  improved  by  regularly  updating  altitude  and 
position.  The  remaining  errors  will,  however,  be  of  a size  which 
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will  not  allow  a useful  geodetic  application  of  the  filtered 
data.  These  findings  from  a probabilistic  error  analysis  were 
confirmed  by  results  obtained  by  Szabo  and  Anthony  (1971)  in 
an  analysis  of  actual  measurements. 

In  the  second  method  structural  differences  between  the 
gravitational  and  inertial  fields  are  used  to  separate  the 
two  effects.  These  differences  show  up  in  the  second  and  higher 
order  gradients  of  the  force  fields.  Therefore,  additional 
measurements  are  necessary  in  this  case.  Moritz  (1967)  has  shown 
that  for  an  aircraft  with  inertial  stabilization  the  second 
derivatives  of  the  force  field  do  not  contain  inertial  distur- 
bances, so  that  purely  gravitational  second-order  gradients  can 
be  measured.  They  are  used  to  obtain  the  gravitational  force 
vector  by  integrating  along  the  flight  path.  It  should  be  noted 
that,  in  contrast  to  the  first  method,  a rigorous  separation  of 
gravitation  and  inertia  is  possible  in  this  case,  and  that  from 
a theoretical  point  of  view  this  approach  is  preferable.  The 
practical  difficulties  originate  in  the  design  of  instruments 
accurate  enough  to  make  an  application  feasible.  Advances  in 
instrument  development  have  been  rapid  during  the  last  years 
and  a gradiometer  with  an  accuracy  of  a few  Eotvos  may  be 
available  in  the  near  future.  Therefore,  the  capabilities  of 
an  airborne  gradiometer  system  are  studied  in  this  report. 

The  accuracy  study  will  be  performed  using  the  method  of 
least-squares  collocation.  There  are  three  reasons  why  this 
approach  seems  to  be  especially  suited  for  the  problem  at  hand. 
First,  it  allows  the  combination  of  heterogeneous  data  in  a 
consistent  way.  This  is  very  important  because  geoidal  heights, 
gravity  anomalies,  and  different  second-order  gradients  will  be 
used  as  measurements.  They  must  be  evaluated  in  such  a way  that 
their  common  origin  from  the  same  anomalous  gravity  field  is 
part  of  the  system.  In  least-squares  collocation  this  is 
achieved  by  describing  the  statistical  structure  of  the  field 
by  a covariance  function.  Second,  mean  gravity  values  at  ground 
level  must  be  estimated  using  point  values  on  a profile  in 


flying  altitude  and  additional  information  on  ground.  This 
involves  interpolation  between  profiles,  downward  continuation, 
combination  of  different  quantities,  and  estimation  of  mean  values. 
All  these  steps  can  be  united  in  a single  step  procedure  in  the 
collocation  method.  This  is  impossible  when  using  the  corres- 
ponding integral  formulas.  Third,  different  assumptions  on  the 
structure  of  the  gravity  field  and  on  the  accuracy  of  the 
measurements  must  be  investigated.  Again  this  is  very  simple 
with  the  collocation  method  because  it  only  involves  a change 
of  the  fundamental  covariance  function  or  of  the  error  variances. 

Basically,  this  is  a deterministic  approach.  Its  results 
are  valid  if  the  assumptions  on  the  structure  of  the  gravity 
field  and  on  the  accuracy  of  the  measurements  are  not  too  far 
from  reality.  The  covariance  functions  have  been  varied  in  a 
rather  wide  range  so  that  they  will  hopefully  satisfy  these 
requ i rements . Therefore, no  investigations  are  made  as  to  the 
influence  of  an  incorrect  covariance  function.  If  necessary, 
results  published  in  Moritz  (1976)  can  be  used  as  a guideline 
for  the  present  case.  The  influence  of  measuring  errors  will 
be  analysed  in  a special  section.  A probabilistic  error  analysis 
has  not  been  attempted  because  adequate  information  on  the  error 
processes  affecting  gradiometer  measurements  is  not  available. 

A short  survey  of  the  following  sections  will  conclude 
this  part.  An  introduction  to  the  mathematical  structure  of  the 
problem  is  given  in  section  2.  Section  3 discusses  the  deter- 
mination of  a spatial  covariance  function  from  empirical  data. 

Since  this  section  is  fundamental  for  the  rest  of  the  report,  the 
arguments  are  presented  in  some  detail.  The  mechanism  of  spatial 
interpolation  using  heterogeneous  data  is  studied  in  sections 
4 and  5.  In  section  6 the  effect  of  measuring  errors  is  dis- 
cussed and  results  are  applied  in  sections  7 and  8.  Questions 
related  to  the  design  of  experiments  are  treated  in  section  7 
as  e.g.  the  spacing  of  flight  profiles  to  obtain  mean  values 
of  certain  accuracies  for  blocks  of  given  sizes.  The  contri- 
bution of  highly  accurate  satellite  altimetry  to  a combined 
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accel erometer-grad i ometer  system  is  considered  in  section  8. 
Finally,  some  numerical  problems  are  summarized  in  section  9. 


2.  Mathematical  Background  and  Numerical  Procedure. 


The  objective  of  this  section  is  to  present  the  main 
formulas  and  to  discuss  their  implications.  This  may  serve  as 
a simplified  introduction  into  the  mathematical  structure  of 
the  problem.  No  derivations  are  given.  They  can  be  found  in 
Moritz  (1971,  1975).  For  a deeper  understanding  Moritz  (1967) 
should  be  consulted. 

It  has  been  indicated  in  the  introduction  why  it  is  not 
possible  to  separate  gravitation  and  inertia  by  using  accelero- 
meter measurements  only.  Let  us  now  consider  in  which  way  the 
use  of  gradiometer  measurements  will  change  the  situation.  This 
may  best  be  seen  from  the  system  of  ordinary  differential  equa- 
tions given  in  Moritz  (1975) 

Tf"1,T.+T('1)T.-T.+f.=0  2.1 

13  j 13  3 i i 

where  T ...  anomalous  gravitational  potential; 

T.  ...  = |^—  i = 1,2,3  first-order  gradients  of  T ; 

i o x . 

3 2 T 

T = — i , j = 1,2,3  second-order  gradients 

1 ] o X . d X . 

i 3 
of  T ; 

f ...  component  of  accelerometer  output  (gravitational 
+ inertial  force)  . 

All  quantities  in  equation  (2.1)  are  regarded  as  functions  of 
time.  The  dots  therefore  denote  differentiation  with  respect 
to  time.  The  summation  convention  has  been  used,  i.e.  an  in- 
dex occuring  twice  in  a product  implies  summation. 

Equation  (2.1)  is  a system  of  linear  second-order  differ- 
ential equations  for  the  gravity  disturbance  vector  T^  . If 


- ..  . *-i-  x 
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the  quantities  and  T i are  given  by  measurement  and  the 

matrix  with  elements  T is  assumed  to  be  invertible,  equa- 
ls 

tion  (2.1)  may  be  solved  by  the  usual  numerical  methods.  Thus, 
it  is  indeed  possible  to  separate  gravitation  and  inertia,  ob- 
taining as  a result  the  purely  gravitational  quantities  T , 
i.e.  the  deflections  of  the  vertical  and  the  gravity  disturbance 
which  usually  is  not  very  different  from  the  gravity  anomaly. 

The  practical  solution  of  equation  (2.1)  may  pose  some 

problems.  First  of  all,  we  need  initial  values  T°  at  time  t 

1 o 

which  may  be  difficult  to  get.  Second,  besides  the  measured 
gradients  T ^ , we  need  their  derivatives  with  respect  to 
time  T.j  . This  involves  a knowledge  of  the  velocity  of  the 
aircraft  which  is  only  approximately  known  because  the  measured 
values  contain  the  influence  of  the  anomalous  gravity  field. 
Therefore,  an  iterative  solution  will  be  advantageous  in  prac- 
tical computations.  The  formulas  for  such  a solution  are  also 
given  in  Moritz  (1975) 

vt(t)  = u ° - / fi(s)ds 
o 

Ti(t)  = T°  + / Tij(s)v.(s)ds 

° t s 

Ti(t)  = T (t)  + / / T (s)T  (r)drds 

s = t r = t J J 

o o 

where  ui  are  the  three  velocity  components  and  the  other 
notations  are  as  defined  in  equation  (2.1).  Again  all  quanti- 
ties are  regarded  as  functions  of  time,  i.e.  the  variables 
t , r , and  s all  refer  to  time  and  are  only  used  to  dis- 
tinguish the  integration  variables  unambigously  from  one 
another.  The  upper  case  ° denotes  the  value  of  the  function 
at  some  initial  time  t 

o 

The  procedure  can  now  be  read  from  the  formulas.  With 
initial  values  u°  and  T°  given,  we  first  compute  the 
approximate  velocity  components  vL  at  time  t . These  values 


2.2 

2.3 

2.4 
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differ  from  the  exact  velocity  components  ujt)  by  the  effect 
of  the  a noma lous  gravity  field  in  the  accelerometer  output  f . 
Assuming  for  a moment  the  gravitational  effect  to  be  absent, 
formula  (2.3)  would  give  us  the  exact  result,  i.e.  the  components 
of  the  gravity  disturbance  vector.  Since  this  is  not  the  case, 
we  get  the  approximation  T ( t ) . Using  this  approximation  in  the 
integral  term  of  equation  (2.4)  will  give  more  accurate  values 
T^t)  which  again  can  be  used  in  the  same  way.  Convergence  pre- 
sumed, we  will  finally  obtain  the  exact  gravity  disturbance  vector. 
It  should  be  noted  that  at  the  same  time  the  exact  velocity 
components  u^t)  can  be  computed.  The  integral  term  in  equa- 
tion (2.4)  has  been  called  interaction  term  because  it  expresses 
the  interaction  between  gravitation  and  inertia.  The  size  of 
this  term  compared  to  T is  important  for  the  speed  of  con- 
vergence of  the  iterative  procedure. 

The  integrations  performed  in  equations  (2.2)  to  (2.4)  are 
all  along  the  flight  path.  From  practical  considerations  this  is 
very  advantageous  because  measurements  can  be  restricted  to  a 
specified  region  and  computations  can  be  performed  in  real 
time.  Difficulties  as  in  Stokes'  or  Veni ng-Mei nesz ' integral 
formulas  where  global  data  coverage  is  required  do  not  occur. 

A simple  analogy  from  classical  geodesy  is  the  a strogeodet i c 
determination  of  the  geoid  where  the  deflections  of  the  vertical 
are  also  needed  along  a profile  only,  in  order  to  compute  the 
geoid  along  this  profile.  It  is  interesting,  however,  that  an 
analogue  to  Stokes  formula  also  exists  for  the  radial  second- 
order  gradient  T^^  . It  has  been  derived  in  Moritz  (1967) 


a 


where  R is  the  mean  radius  of  the  earth,  tp  is  the  spherical 
distance  and  S^ <p ) is  an  analogue  to  Stokes  function  and 
given  by 
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Sx(*)  = (2-6sin2J) 1 n(l+l/sin|)  -3+6  s i n|  . 2.6 

The  integration  in  (2.5)  has  to  be  extended  over  the  surface 
of  the  earth.  A general  use  of  this  formula  cannot  be  re- 
commended because  the  effect  of  the  remote  zones  is  even 
stronger  than  in  Stokes'  formula.  But  if  a good  coverage  of 
T^-values  is  given  in  a certain  area  and  enough  outside  in- 
formation is  available  from  other  sources,  a method  similar 
to  that  of  Marsh  and  Vincent  (1974)  may  be  applied  to  compute 
the  fine  structure  of  the  geoid  in  that  area.  Such  a procedure 
would  also  allow  interesting  comparisons  with  the  collocation 
procedure  to  be  treated  next. 

So  far,  we  have  discussed  the  problem  of  separating 
gravitation  and  inertia  without  reference  to  the  quantities  we 
actually  want  to  obtain  by  our  computations.  Usually,  these 
quantities  will  be  mean  gravity  values  in  blocks  of  specified 
sizes  on  the  surface  of  the  earth.  To  derive  them  from  the 
measured  or  the  integrated  values,  downward  continuation  and  avera- 
ging procedures  must  be  applied.  This  should  be  done  in  such  a 
way  that  no  relevant  information  of  the  measuring  process  is 
lost  in  the  subsequent  computations.  This  is  difficult,  however, 
with  classical  procedures.  Taking  as  an  example  formulas  (2.2) 
to  (2.4),  we  will  obtain  the  three  components  of  the  gravity 

disturbance  vector  T , T , T . Downward  continuation  of 

0 X r 

T , as  that  of  Ag  , does  not  involve  T or  T,  . Thus, 

r <j>  A 

valuable  information  is  not  used  at  all.  Similarly,  downward 

continuation  of  the  T -gradient  would  not  make  use  of  the 

r r 

four  other  independent  second-order  gradients.  The  situation 
becomes  even  more  complicated  if  measurements  from  other  sources 
which  are  of  different  type  must  be  combined  with  the  gradio- 
meter  data.  Therefore,  we  need  a method  which  permits  the  combi- 
nation of  heterogenous  data  in  different  altitudes  and  their 
downward  continuation  to  ground  level.  To  obtain  mean  values 


from  point  values  an  interpolation  procedure  must  be  available 
which  uses  information  about  the  average  behaviour  of  the  quan- 
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t i t i es  involved. 

Moritz  (1971)  has  proposed  to  use  1 ea s t- squa res  collo- 
cation for  this  purpose.  The  computational  formulas  are  very 
simple,  namely 

s = C C” 1 x 2.7 

SX  XX 

and 

E = C - C C"1CT  2.8 

ss  ss  sxxxsx 

where  s , the  signal,  is  the  quantity  to  be  estimated,  e.g. 

gravity  anomalies;  x is  the  vector  of  observations,  e.g. 

second-order  gradients,  deflections  of  the  vertical,  gravity 

anomalies,  height  anomalies;  Cgx  is  the  covariance  matrix 

connecting  signal  and  observations;  Cxx  is  the  autocovariance 

matrix  of  the  observations;  C is  the  autocovariance  matrix 

s s 

of  the  signal.  Given  a number  of  observed  values  and  an  appro- 
priate covariance  function,  the  signal  quantities  s and  their 
error  covariance  matrix  E may  be  estimated  from  formulas 
(2.7)  and  (2.8).  Since  we  are  interested  in  accuracy  estimates 
only  formula  (2.8)  will  be  used  in  the  sequel. 

The  concept  of  a covariance  function  plays  a decisive 
role  in  the  collocation  method.  Therefore  it  will  be  discussed 
in  detail  in  section  3 and  we  can  restrict  ourselves  to  some 
more  general  remarks  at  this  point.  Basically,  the  covariance 
function  supplies  information  on  the  structure  of  the  gravity 
field.  As  all  measured  quantities  belong  to  the  same  gravity 
field  and  are  related  by  functional  equations,  it  is  possible, 
under  certain  assumptions,  to  derive  their  different  covariance 
functions  from  one  basic  function.  This  is  the  reason  why  the 
estimation  by  formuals  (2.7)  and  (2.8)  can  combine  heterogenous 
data  in  a consistent  way.  Furthermore,  the  basic  covariance 
function  can  be  selected  in  such  a way  that  it  can  be  analy- 
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tically  continued  down  to  sea  level  and  all  difficulties 
usually  encountered  in  downward  continuation  can  be  avoided. 
Interpolation  and  downward  continuation  are  then  solved  by 
the  same  algorithm.  The  information  on  the  different  altitudes 
is  contained  in  the  covariance  function  since  tie  height 
difference  is  one  variable  of  this  function.  The  estimation  of 
mean  values  is  also  possible  by  means  of  the  covariance  function. 
The  averaging  process  is  applied  to  the  covariance  function  itself, 
changing  the  point  anomaly  function  to  a mean  anomaly  function, 
and  again  the  algorithm  of  formulas  (2.7)  and  (2.8)  is  appli- 
cable. If  measuring  errors  must  be  taken  into  account  this  can 
be  achieved  by  adding  their  effect  to  the  autocovariance  matrix 
of  the  measurements.  Thus,  all  the  different  steps  can  be  per- 
formed by  appropriate  changes  of  the  covariance  function,  making 
the  estimation  of  mean  values  from  point  values  of  heterogeneous 
data  in  different  altitudes  a single  step  procedure.  This  makes 
the  method  so  attractive  for  applications. 

Besides  the  signal  which  can  be  regarded  as  a stochastic 
quantity  systematic  effects  may  be  taken  into  account.  The 
corresponding  formulas  are 


X = (AtC‘1A)'1AtC  1 x 

' XX  7 XX 


s = C C A(x  - AX) 

SX  XX  ' 9 


Exx  = (aV^A)-1 


E = C - C C 1CT  + C C-1AE  ATC“1CT 

SS  SS  SX  XX  SX  SX  XX  XX  XX  SX 


where  X are  the  systematic  parameters  to  be  estimated,  A 
is  the  coefficient  matrix  of  the  parameters,  and  Exx  is 
their  error  covariance  matrix.  Simplifications  of  these 
formulas  are  possible  in  special  cases,  see  Schwarz  (1974). 
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3^ Choice  of  a Covariance  Function. 


It  has  been  pointed  out  that  in  least-squares  collocation 
the  information  on  the  structure  of  the  gravity  field  is  sup- 
plied by  the  covariance  function.  The  determi nat i on  of  such  a 
function  from  empirical  data  is  usually  governed  by  two  view- 
points. On  the  one  hand  the  covariance  function  should  re- 
present certain  statistical  properties  of  the  gravity  field  as 
determined  from  the  data;  on  the  other  hand  it  should  be  of  a 
simple  analytical  form  so  that  the  necessary  derivations  can 
be  performed  in  such  a way  as  to  secure  computational  efficiency. 
The  second  requirement  has  so  far  only  been  met  by  isotropic 
and  homogeneous  covariance  functions,  i.e.  functions  which  are 
invariant  with  respect  to  rotations  and  transl at  ions . There- 
fore the  following  discussion  will  be  restricted  to  these 
functions. 

A question  which  arises  naturally  when  fitting  empirical 
data  to  some  kind  of  model  function  is:  how  many  parameters  are 
essential  for  the  model  function?  This  question  has  recently 
been  studied  by  Moritz  (1976)  for  isotropic  and  homogeneous 
covariance  functions  in  the  plane  which  can  be  extended  into 
outer  space.  He  has  shown  that  rotationally  symmetric  harmonic 
covariance  functions  in  space  usually  have  planar  equivalents  and 
that  results  obtained  for  the  planar  approximations  carry  over 
to  the  spherical  case  with  only  small  modifications.  These  re- 
sults will  be  used  in  the  sequel. 

In  Moritz  (1976)  three  essential  parameters  are  given  for 
a covariance  function  C(s)  : the  variance  , the  corre- 
lation length  £ , and  the  curvature  parameter  x • Their 
geometrical  interpretation  is  given  by  fig.  3.1.  The  variance 

C is  the  value  of  the  covariance  function  C(s)  for  the 

o 

argument  s = 0 

Cq  = C ( 0 ) . 3.1 


« 


C(si 


Covariance 

function 


Fig.  3.1  Geometrical  i nt er preta t i on 
of  essential  parameters 


The  correlation  length  5 is  the  value  of  the  argument 
for  which 


The  curvature  parameter  x is  related  to  the  curvature 
of  the  covariance  curve  at  s = 0 by 


where  k is  the  curvature.  If  we  norm  the  covariance  curve 
by  setting  = 1 and  5=1  we  obtain 


which  explains  the  name  chosen. 


C ( 5 ) = C o / 2 . 
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The  three  parameters  have  been  called  essential.  This 
does  not  mean  that  two  covariance  functions  which  have  the 
same  numerical  values  for  the  parameters  are  equal.  But  their 
general  behaviour  with  respect  to  interpolation  will  be  very 
similar,  i.e.  the  interpolation  errors  will  almost  be  the  same 
in  both  cases.  The  variance  C is  a kind  of  scale  factor  for 

o 

the  interpolation  errors.  Two  covariance  functions  which  differ 
only  in  their  variances  C;  and  C;  will  have  interpolation 
errors  and  m^  related  by 

= r~2  m2  • 3-5 

Relations  between  the  other  parameters  cannot  be  put  into  such 
a simple  form  but  it  can  generally  be  said  that  x character- 
izes the  behaviour  of  the  covariance  function  for  small  dis- 
tances s , while  £ describes  the  behaviour  for  distances 
on  the  order  of  £ itself.  This  indicates  that  the  shape  of 
the  covariance  curve  at  distances  larger  than  about  1.5  ? 
is  not  very  decisive,  a fact  which  is  confirmed  by  numerical 
computations.  On  the  other  hand  a realistic  choice  of  x is 
very  important  because  the  curvature  at  the  origin  influences 
the  interpolation  for  small  distances.  A model  which  is  only 
fitted  to  Cq  and  C ( ^ ) may  produce  a strongly  distorted 
picture  of  the  actual  covariance  behaviour.  It  can  e.g.  be 
shown  that  for  a reasonable  overall  fit  of  the  empirical  gravity 
anomaly  covariance  function,  the  variance  G of  the  horizontal 
gradients  of  Ag  may  differ  by  a ratio  of  1:30  depending  on 
the  model  chosen.  Since  this  variance  is  related  to  the  curvature 
parameter  x by 

G = x • C H2  , 3.6 

o o 

see  (Moritz,  1976  , p . 2 5 ) , the  necessity  of  a third  essential 
parameter  is  quite  evident.  In  practice,  the  curvature  at  the 
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the  origin  is  usually  not  very  well  known  from  empirical  data. 

On  the  other  hand,  the  variance  G is  more  accessible  to 

o 

statistical  estimation.  It  is  therefore  advisable  to  select 
the  basic  covariance  model  using  G as  the  third  parameter. 

The  curvature  parameter  \ may  then  be  determined  by 

X = G • e2/C  . 3.7 

O O 

Thus,  we  will  regard  C , C > and  G as  the  three  essential 

o o 

parameters  in  the  sequel. 

As  can  be  seen  from  this  example  the  choice  of  an  ade- 
quate covariance  model  is  influenced  by  the  amount  and  type  of 
empirical  data  available.  Presently,  good  statistical  in- 
formation on  some  important  parameters  is  1 ac k i ng  . Furt hermore , 
none  of  the  covariance  models  proposed  so  far  will  fit  all 
empirical  data  equally  well.  Probably  this  can  only  be 
achieved  by  combining  different  models.  It  seems  therefore 
reasonable  to  distinguish  global  and  local  applications  at 
this  point.  It  must  be  stressed,  however,  that  this  distinction 
is  purely  from  practical  reasons.  In  global  applications  a good 
fit  to  the  covariance  function  of  mean  gravity  anomalies  of 
certain  block  sizes  and  to  the  low  degree-variances  as  deter- 
mined from  satellites  is  satisfactory.  In  local  applications 
a good  fit  to  a point  anomaly  function  and  to  the  variance  of 
its  gradients  is  necessary.  It  is  therefore  not  advisable  to 
derive  local  covariance  functions  from  a global  function  deter- 
mined in  the  above  way.  Mean  gravity  anomalies  do  not  contain 
reliable  information  on  the  curvature  parameter  of  the  corres- 
ponding point  covariance  function.  If  the  variances  Cq  or 
G^  are  fitted  to  their  empirical  values  afterwards  the  in- 
correct information  is  already  part  of  the  model  and  will  in- 
fluence all  subsequent  computations.  Consequently,  the  para- 
meters essential  for  a local  covariance  function  must  already 
be  considered  when  selecting  the  model. 

H 


In  the  following  only  local  covariance  functions  are 
discussed  because  the  use  of  grad i ometer  data  makes  a realistic 
local  behaviour  much  more  important  than  a good  global  fit. 
Therefore,  we  need  as  essential  parameters  the  variance  Cq 
of  the  point  gravity  anomalies,  the  correlation  length  £ 
of  the  correspondi ng  covariance  function  and  the  variance  G 
of  the  horizontal  gradients  of  the  gravity  anomalies.  Further- 
more, since  aerial  gradiometry  has  important  advantages  when 
applied  over  sea,  differences  in  the  gradient  behaviour  over 
sea  and  over  land  should  be  considered.  As  a typical  example 
we  will  regard  a sea  surface  area  with  depth  between  3 and 
6 km.  This  depth  range  covers  more  than  50  % of  the  surface  of 
the  earth,  i.e.  more  than  70  % of  the  total  sea  surface. 

Representative  covariance  functions  of  point  gravity 
anomalies  have  not  been  published  in  the  last  decade;  the 
careful  investigations  of  Tscherning  and  Rapp  (1974)  have  been 

O O ^ 

devoted  to  a covariance  function  from  1 x 1 mean  gravity 
anomalies,  and  the  point  anomaly  function  derived  from  it  will 
not  be  used  here  for  reasons  mentioned  above.  Therefore,  the 
basic  material  to  obtain  estimates  for  C and  £ has  been 

O 

taken  from  publications  of  Hirvonen  (1962)  and  Kaula  (1966). 
Hirvonen's  covariance  function  is  a purely  local  one  deter- 
mined from  12  samples  (area  1°  x 1°)  in  the  state  of  Ohio. 
Kaula's  covariance  function  is  more  representative  coming 
from  9 local  (area  2°  x 2°)  and  8 regional  samples  (area 
10°  x 10°),  half  of  them  ocean  areas.  Using  these  empirical 
data,  two  local  covariance  functions  will  be  derived  which 
should  be  different  enough  to  realistically  cover  the  range 
of  possible  functions.  The  information  needed  in  the  sequel  is 
summarized  in  table  3.1: 
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A reliable  estimate  of  the  variance  G is  difficult  to  get. 

° 2 

The  values  found  in  different  publications  vary  between  30  E 
and  330  EJ  for  the  horizontal  gradients  of  fig  , where 

1 E = 10~9sec”2  = 0.1  mgal/km  . 

Most  of  these  values  are  not  very  representative  because 
sample  sizes  are  small  and  are  taken  from  one  area  only. 
Furthermore,  all  of  them  are  from  terrestrial  measurements  and 
are  influenced  by  topographical  effects.  Sea  surface  gradients 
for  the  given  depth  range  will  probably  have  a much  smoother 
behaviour  than  their  terrestrial  equivalents.  There  is  another 
reason  why  for  aerial  gradiometry  the  parameter  G may  be 
expected  close  to  the  lower  limit  of  the  above  estimates.  The 
measurements  are  not  actually  point  values  but  mean  values  at 
interval!  of  10  sec,  i.e.  at  a distance  of  about  2.6  km  for  a 
flying  speed  of  500  knots.  Naturally,  these  mean  values  will 
be  smoother  than  point  values.  Assuming  therefore 

30  E2  - G - 200  E2  , 3.8 

O 

we  obtain  from  table  3.1  and  formula  (3.3) 

1.4  - Xj  - 8.5  1.7  - x2  - 10.8  3.9 

where  the  indices  1 and  2 refer  to  the  Hirvonen  and  the  Kaula 
estimates  respectively. 

With  estimates  for  the  three  essential  parameters  given. 


we  can  select  a covariance  model  which  will  meet  the  require- 
ments. A rotational ly  symmetric  harmonic  covariance  function 
in  space  can  be  expressed  in  the  general  form 


where  P and  Q are  two  points  in  space  with  radius  vectors 
rp  and  r^  respectively,  ^ is  the  angle  between  rp  and 
r . , R is  the  radius  of  some  sphere  - frequently  called  the 
Bjerhammar  sphere  - completely  inside  the  mean  terrestrial 
sphere  (R  = 6371  km),  P (cosj')  are  the  Legendre  polynomials, 
and  k are  positive  coefficients.  We  will  use  the  substi- 

u 

tu  t i on 


s 


r r 
p Q 


3.11 


since  no  confusion  with  the  distance  variable  s can  occur  in 
the  sequel  and  we  obtain 

oo 

K(P,Q)  = A l kns“  + 1Pii(cos*)  . 3.12 

11=0 

Different  models  of  the  covariance  function  can  be  obtained 
by  defining  k in  different  ways.  Some  of  them  are  given 
in  table  3.2: 


Name 

k 

n 

X of  planar 
equivalent 

Reciprocal  Distance  C.F. 

i 

3.00 

Poi sson  C.F. 

2n  + 1 

1.76 

Logarithmic  C.F. .type  1 

1/n 

variable 

" " .type  2 

1/  (n-1)  (rt-2) 

II 

" .type  3 

1 / ( n- 1 ) ( n- 2 ) ( n + B ) 

II 

Table:  3.2:  Covariance  function  models  according  to  equation  (3.11) 


The  names  have  been  taken  from  Moritz  (1976)  where  also  the 
determination  of  the  curvature  parameter  x for  each  function 
is  discussed  in  detail.  The  logarithmic  covariance  models  of 
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type  2 and  type  3 have  been  used  by  Tscherning  and  Rapp  (1974) 
to  obtain  a global  covariance  function  from  empirical  data. 

The  reciprocal  distance  covariance  function  and  the  Poisson 
covariance  function  have  already  been  considered  by  Krarup 
(1969).  Since  the  information  on  the  empirical  x“Va^ues  is 
only  given  in  a certain  range  it  seems  advisable  not  to  use  the 
first  two  models  because  we  may  force  the  data  by  keeping  a 
fixed  value  for  x • On  The  other  hand,  we  may  arrive  at  un- 
realistically large  x'values  when  using  one  of  the  logarithmic 
covariance  functions  (see  Moritz,  1976,  p.48). 

Let  us  therefore  first  consider  the  parameters  in  equa- 
tion (3.12)  and  their  relation  to  the  three  essential  para- 
meters Cq  , 5 , and  x(Gq)  • The  value  of  K(P,Q)  can  be 
varied  in  three  ways  once  the  model  for  the  coefficients 

has  been  selected.  We  can  change  the  values  of  A and  s 
on  the  right-hand  side  of  equation  (3.12),  and  we  can  substract 
a certain  number  of  coefficients  k . A change  of  A will 

II 

only  change  the  scale  of  the  covariance  function,  i.e.  A is 
directly  related  to  C . If  s becomes  sma ller,  the  covar- 

O 

iance  curve  flattens  out,  i.e.  e.  becomes  larger  and  Cq 
becomes  smaller.  Since  G decreases  even  more  rapidly  than 

C , the  parameter  x stays  almost  constant  for  the  loga- 

o 

rithmic  covariance  function  of  type  2 and  decreases  for  type  3. 

If  we  subtract  the  first  m coefficients  k the  covariance 

n 

curve  becomes  steeper  and  the  parameters  c,  , Cq  and  x 
become  smaller.  The  value  of  G^  does  not  change  very  much 
as  long  as  we  stay  in  the  low  frequency  range.  This  shows  that 
the  local  character  of  the  second-order  gradients  is  well 
represented  by  these  models.  In  the  practical  determi na t i on 
it  is  most  important  to  determine  first  an  appropriate  range 
of  values  for  s , because  we  may  otherwise  obtain  unrealistic 
values  for  G and  x • This  can  best  be  done  by  fixing  the 

ratio  C : G because  it  is  most  sensitive  to  changes  in  s . 

o o 

Thereafter,  j;  can  be  fitted  by  subtracting  a certain  number 
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of  coefficients  and  finally  the  function  is  scaled  by  selec- 
> ting  an  appropriate  value  for  A . Table  3.3  shows  the  inter- 

play between  C , G » C » x and  s . 


Logarithmic  cov. 

f ct . , 

Logarithmic  cov. 

f ct . , 

type  2 

type  3 

s 

C : G 

0 0 

C 

X 

C : G 
0 0 

5 

X 

0.999 

0.4:1 

11.2 

3.02 

2.7:1 

18.12 

0.998 

1.6:1 

22.2 

3.03 

9.0:1 

11.93 

0.997 

3.6:1 

33.4 

3.04 

17.8:1 

9.49 

0.996 

6.5:1 

44.6 

3.05 

28.8:1 

153.8 

8.21 

0.995 

10.2:1 

55.8 

3.06 

41.8:1 

174.7 

7.31 

0.994 

14.7:1 

67  . 1 

3.07 

56.3:1 

197.6 

6.93 

0.993 

19.9:1 

78.3 

3.07 

72.5:1 

216.3 

6.45 

0.992 

26.0:1 

89.5 

3.08 

90.2:1 

234.6 

6.10 

Table  3.3  Covariance  models  and  essential  parameters. 

The  best  estimate  of  a global  point  gravity  anomaly 
variance  is  the  value  C = 1795  mgal2  determined  by  Tscherning 

O 

and  Rapp  (1974).  Using  table  4 of  that  publication,  we  get  a 

2 

somewhat  lower  C = 1576  mgal  for  ocean  areas  in  the  3 to 

o 

6 km  depth  range.  To  get  an  estimate  comparable  to  Kaula's 
sampling  areas  we  have  to  subtract  the  influence  of  the  low 

frequencies  with  half  range  larger  than  10°.  Using  satellite 

12 

deri.ed  coefficients,  we  get  values  of  about  1590  mgal  and 

1375  mgal  for  the  above  estimates.  Thus,  Cq  = 1500  mgal  has 
been  chosen  as  a round  value  for  a covariance  function  of  Kaula 
type.  The  effect  of  a change  in  the  variance  Cq  on  the  mean 
square  error  can  then  easily  be  computed  by  formula  (3.5). 
Considering  the  inequality  (3.8),  the  ratio  C^  : G should 
be  about  10:1  to  15:1.  Using  this  ratio  to  fix  the  range  of  s , 
the  correlation  length  5 and  the  scale  of  the  covariance 
function  have  been  determined  in  the  way  described  above. 
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.. 


- - — ~ 


o empirical  values 
— covariance  model 


Fig. 3. 2 Fit  to  Hirvonen's  empirical  covariances 


• empirical  values 
— covariance  model 


Fig. 3. 3 Fit  to  Kaula's  empirical  covariances 
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Similarly,  Hirvonen's  data  have  been  fitted  using  a C -value 

2 ^ 

of  500  mgal  in  order  to  avoid  overly  optimistic  interpolation 

errors.  The  results  of  the  curve  fitting  is  shown  in  figures 

3.2  and  3.3  with  C -values  normed  to  C = 1 . In  both  cases 

o o 

a logarithmic  covariance  function  of  type  2 has  been  used.  A 
function  of  type  3 may  give  a better  agreement  for  the  low 
satellite  derived  degree  variances  but  for  a local  covariance 
function  this  is  not  very  important.  On  the  other  hand,  a 
function  of  type  2 is  somewhat  simpler  to  handle  numerically. 
There  is  a good  overall  fit  for  both  curves  up  to  distances  of 
2 £ • As  has  already  been  pointed  out  deviations  of  the  curves 
for  larger  distances  do  not  affect  the  results  of  interpolation 
A number  of  characteristic  parameters  for  both  functions  are 
given  in  table  3.4  where  the  Hirvonen  type  function  has  been 
called  "local"  and  the  Kaula  type  function  "regional". 


Parameter 


Variance  of  gravity 
a nomal i es  C 


Variance  G 

O 

of  the  horizontal  gra- 
dients of  a ~ 


Local  Regional 

covariance  function  covariance  function 


500  mgal 


52.5  E 


1500  mgal 


111.4  E 


9.5  : 1 

13.5  : 1 

Correlation  length  £ 

43  km 

61  km 

Curvature  parameter  x 

1 . 9 

2.8 

s = 

0.994 

0.994 

Number  of  coefficients 

70 

12 

subtracted 

Table  3.4  Characteristic  parameters  of  the  two  covariance  func 
tions  selected. 
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Using  the  above  parameters  the  covariance  functions 
of  the  geoidal  undulations  K(u),  of  the  gravity  anomalies 
C(ip)  .and  of  the  second-order  radial  gradients  G ( ip ) have 
been  computed  for  altitudes  of  0 km  and  10  km.  They  are 
shown  in  figures  3.4  and  3.5.  It  should  be  noted  that  the 
scale  on  the  ^ - a x i s is  different  for  the  three  cases  but  equal 
for  the  corresponding  local  and  regional  covariance  functions. 
There  are  typical  similarities  and  distinctions  between  the 
two  types  of  functions.  Similar  features  are  produced  by  the 
order  of  differentiation  which  in  fig.  3.3  and  3.4  increases 
from  top  to  bottom.  The  correlation  length  £ reduces  with 
each  differentiation  while  the  attenuation  with  growing  alti- 
tude becomes  stronger  with  each  step.  Differences  are  also 
quite  evident.  The  more  local  character  of  the  first  covariance 
function  is  apparent  from  the  smaller  £-values  for  K ( ) and 
C(;)  and  from  the  smaller  variances  of  all  three  functions. 

Also  the  reduction  of  £ with  each  differentiation  is  quite 
different.  While  for  the  regional  curve  the  correlation  length 
reduces  by  a factor  of  almost  5 with  each  step,  the  effect  is 
much  smaller  for  the  local  function.  Furthermore,  the  small 
variance  and  the  small  correlation  length  of  the  local  K ( ip ) - 
function  indicate  that  the  variation  of  the  geo  id  is  considered 
in  a very  limited  area  only. 

These  features  can  be  used  to  get  a first  impression 
about  the  effect  of  different  types  of  data  on  the  estimation. 
Loosely  speaking  the  correlation  length  £ gives  the  zone  of 
influence  of  the  quantity  described  by  the  covariance  function. 
Thus,  altimeter  data  will  have  a noticeable  effect  on  the 
estimation  up  to  a distance  of  3°. 5,  gravity  data  up  to  about 
0°.6,  and  gradiometer  data  up  to  about  0°.15.  On  the  other  side, 
altimeter  data  will  introduce  much  stronger  correlations  between 
neighbouring  points  than  gravity  or  gradiometer  data.  Further- 
more, the  influence  of  different  types  of  data  will  be  more 
distinct  for  the  regional  than  for  the  local  covariance  function. 
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These  features  can  advantageously  be  used  for  the  combination 
of  heterogeneous  data  and  for  the  design  of  stable  numerical 
procedures,  although  for  a full  assessment  crosscorrelations 
have  to  be  taken  into  account. 

4.  Interpolation  at  Flight  level. 

The  first  tentative  conclusions  drawn  at  the  end  of  the 
preceding  section  will  now  be  substantiated  by  detailed  compu- 
tations. The  basic  formula  has  already  been  given  in  section  2 

E = C - C C-‘CT  . 4.1 

ss  ss  sxxxsx 

The  values  to  be  estimated,  the  components  of  the  signal  s , 
are  .'.g-values  in  arbitrary  positions  between  the  profiles.  E^ 
is  their  error  covariance  matrix  and  the  square  roots  of  the 
diagonal  terms  of  this  matrix  are  the  standard  errors  of  the 
gravity  anomalies.  The  measurements  x used  in  this  section 
are  gradiometer  data  in  a <j>  , \ , r-system  and  Ag-values  ob- 
tained by  integrating  along  the  profile  according  to  formulas 
(2.2)  to  (2.4).  'lean  values  at  intervals  of  10  sec.  are  con- 
sidered as  gradiometer  measurements.  The  elements  of  the  co- 
variance  matrices  are  determined  by  applying  functions  similar 
to  equation  (3.12)  for  two  specific  points  P and  Q or, 
since  we  are  using  isotropic  covariance  functions,  by  intro- 
ducing the  spherical  distance  4)  between  the  two  points  and 
their  respective  heights.  If  measuring  errors  must  be  taken 
into  account,  the  covariance  matrix  of  the  measurements  C 

XX 

has  to  be  augmented  by  the  covariance  matrix  of  the  noise 
C and  we  shal 1 write 

un 

C = C + C 4.2 

xx  mi 
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so  that  formula  (4.1)  reads 

E = C - C C'1CT  . 4.3 

ss  ss  sx  sx 

This  formula  will  always  be  applied  in  the  sequel  for  reasons 
of  numerical  stability  as  discussed  in  section  9.  We  have  used 
standard  errors  of  + 1 mgal  and  + 1 E for  the  gravity  anomalies 
and  the  gradiometer  data  respect i vel y . If  we  regard  these  errors 

as  independent  we  have  a C -matrix  of  the  form 

n ji 

C =o.I  4.4. 

1111  1 

where  o , are  the  error  variances  of  the  different  types  of 
data  and  I is  the  unit  matrix.  Usually  the  differences  between 
this  model  and  pure  interpolation  with  o.  = 0 are  negligible. 

If  necessary , has  been  appropr i ately  reduced. 

It  will  first  be  investigated  how  the  data  can  be  com- 
bined in  an  optimal  way.  There  are  four  obvious  alternatives. 
Following  the  classical  procedure  we  can  regard  Ag  only  and 
interpolate  between  the  profiles.  In  this  case  a large  part  of 
the  existing  information  is  neglected.  Next,  the  use  of  all  three 

first-order  gradients  T.  , T,  , T will  be  considered,  im- 

<J>  a r 

plying  that  the  information  of  five  independent  gradiometer 
measurements  can  be  condensed  into  three  first-order  gradients. 
Alternatively,  we  will  use  five  second-order  gradients  only. 
Finally,  a combination  of  Ag  and  different  second-order  gra- 
dients will  be  used.  It  should  be  noted  that  in  all  cases  where 
integrated  values  are  employed  correlations  via  the  measuring 
errors  are  introduced.  A rigorous  model  should  take  these  corre- 
lations into  account.  Since  their  influence  is  negligible  as  long 
as  the  error  variances  are  small  compared  to  the  signal  variances 
the  simple  error  model  (4.4)  has  been  kept. 

When  considering  formula  (4.1)  it  can  be  seen  that  the 
inversion  of  the  matrix  C is  the  time-consuming  part,  numeri- 
cally. The  dimension  of  C is  given  by  the  number  of  measure- 
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merits  used.  Since  at  each  point  up  to  six  observations  are 
given  and  since  the  points  are  very  dense  along  the  profile, 
the  question  how  to  reduce  the  number  of  measurements  without 
impairing  the  validity  of  the  results  is  not  a trivial  one.  After 
some  preliminary  computations  it  was  found  that  the  arrangement 
shown  in  fig.  4.1  was  satisfactory  to  answer  most  questions  of 
interest.  In  this  case  the  appropriate  range  of  observations 
along  the  profile  should  be  about  three  times  the  distance  from 
the  estimated  point  to  the  profile.  The  arrangement  should  be 
symmetrical  with  respect  to  the  estimated  point.  The  required 
density  of  measurements  on  the  profile  can  be  judged  from  fig. 4. 2 
where  the  number  of  points  in  the  optimal  interval  has  been 
varied.  As  can  be  seen  the  standard  errors  decrease  only  slightly 
with  an  increasing  number  of  points  and  even  an  a rr a ng ement  with 
three  points  on  each  profile  will  bring  reliable  results.  Thus, 
a dense  sequence  of  points  along  the  profile  which  is  necessary 
when  integrating  first-order  gradients,  is  not  needed  for  inter- 
polation purposes.  It  should  be  noted,  however,  that  these  re- 
sults do  not  carry  over  to  the  estimation  of  mean  values.  Further 
more,  to  reduce  the  numerical  work,  only  east-west  profiles  have 
been  considered.  Conclusions  for  the  general  case  can  be  derived 
without  difficulties. 

Results  are  summarized  in  fig.  4.3  for  the  local  covariance 
function  and  in  fig.  4.4  for  the  regional  covariance  function. 

The  standard  error  of  the  interpolated  point  is  plotted  as  a 
function  of  the  profile  spacing.  The  point  is  always  situated 
in  the  middle  between  the  two  profiles.  Asymmetrical  arrange- 
ments have  also  been  investigated  but  are  not  shown  here.  The 
flight  level  is  10  km  above  sea  level.  There  are  two  obvious 
conclusions.  First,  the  accuracy  of  the  interpolated  point  is 
strongly  dependent  on  the  distance  to  the  nearest  profile.  Second 
the  use  of  gravity  anomalies  and  gradiometer  data  will  give  the 
best  interpolation  results.  A more  detailed  examination  of  the 
curves  shows  that  results  are  very  poor  when  using  gradiometer 
data  only.  This  indicates  that  we  can  only  expect  the  fine 


Fig. 4. 3 Local  covariance  function 

Standard  error  of  interpolation 


Fig. 4. 4 Regional  covariance  function 

Standard  error  of  interpolation 
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structure  of  the  gravity  field  from  second-order  gradients.  More 
information  is  needed  to  obtain  good  regional  or  global  values. 
Again,  the  analogy  with  the  astrogeodetic  determination  of  the 
geoid  may  help  to  understand  the  situation.  By  that  method  we 
get  differences  of  the  geoidal  undulations  very  accurately  but 
the  accuracy  of  the  undulations  themselves  will  be  poor  without  a 
well  determined  initial  value.  Similarly,  the  initial  values  in 
equation  (2.3)  play  a decisive  role.  As  soon  as  one  Ag-value  is 
introduced  as  an  observation  results  improve  drastically  even 
if  the  accuracy  of  the  Ag-value  is  not  very  good.  On  the  other 
hand,  when  using  gravity  anomalies  only,  the  detailed  information 
contained  in  the  second-order  gradients  is  missing  and  the  accu- 
racy is  significantly  lower  than  in  the  combination  solution.  The 
same  argument  holds  in  case  of  the  three  first-order  gradients 
for  small  profile  spacings.  For  profile  spacings  larger  than 
1°.5  first-order  gradients  will  give  the  best  solution  but  stan- 
dard errors  are  too  big  to  allow  a useful  geodetic  application 
of  these  large  spacings.  The  best  overall  solution  is  obtained 
by  combining  Ag  and  five  gradients.  But  the  solution  using 
ag  , T , and  T is  almost  as  good.  This  indicates  that  in 
our  special  case  most  of  the  relevant  information  is  contained 
in  these  data.  When  using  north-south  flights  the  corresponding 
quantities  would  be  ag  , T and  T , and  for  flights  with 
an  arbitrary  azimuth  all  four  second-order  gradients  must  be 
used  to  get  a solution  of  comparable  accuracy.  Some  other  data 
combinations  have  been  tried  but  for  the  given  number  of  gra- 
dients the  results  shown  here  are  the  best. 

The  differences  between  the  results  of  the  two  covariance 
functions  are  as  expected.  The  standard  errors  are  comparable 
when  taking  a scaling  factor  according  to  equation  (3.5)  into 
account.  The  more  local  character  of  the  first  function  is 
visible  in  the  fact  that  the  standard  errors  level  off  for  large 
profile  distances.  For  small  profile  distances  the  corresponding 
results  from  both  functions  are  very  close  as  should  be  expected 
from  a good  modeling. 
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5.  Downward  Continuation. 

Downward  continuation  of  gravity  is  an  improperly  posed 
problem  and  the  numerical  difficulties  which  arise  when  looking 
for  an  approximate  solution  have  been  described  in  Schwarz  (1973). 
Downward  continuation  of  second-order  gradients  is  an  even  more 
sensitive  process  because  the  attenuation  of  the  values  with 
height  is  much  stronger.  For  a detailed  discussion  of  the  spectral 
properties  of  such  a process  reference  is  made  to  Rummel  (1975). 
The  mathematical  instability  of  the  problem  is  eliminated  in  the 
collocation  procedure  by  selecting  a covariance  function  that 
can  be  analytically  continued  down  to  sea  level.  In  this  way 
downward  continuation  can  be  regarded  as  a problem  of  spatial 
interpolation.  The  numerical  instability,  i.e.  the  possibility 
that  the  resulting  covariance  matrix  is  close  to  an  ill-con- 
ditioned one,  is  controlled  by  using  a formula  of  type  (4.3) 
where  the  diagonal  of  the  matrix  to  be  inverted  is  augmented  by 
small  amounts.  Thus,  there  is  no  change  in  the  numerical  pro- 
cedure compared  to  the  preceding  section.  The  transition  from 
a plane  to  a spatial  covariance  function  will  only  be  visible  in 
the  covariance  matrices  C and  C because  the  signal  is 

s s sx 

now  at  ground  level  (h  = 0 km)  while  the  measurements  remain 
at  flight  level  (h  = 10  km). 

Again,  some  computations  have  been  made  to  determine  the 
range  and  the  density  of  the  observations.  Results  are  similar 
to  those  of  the  preceding  section  and  no  change  in  the  basic 
design  is  necessary.  As  can  be  seen  from  figures  5.1  and  5.2 
spatial  interpolation  follows  the  same  pattern  as  planar  inter- 
polation if  the  flying  altitude  is  small  compared  to  the  profile 
spacing.  Again,  the  accuracy  of  an  interpolated  point  is  de- 
pendent on  the  distance  to  the  closest  observation  point  and 
the  solutions  combining  Ag-values  and  gradiometer  data  will 
give  the  best  results.  Again,  for  east-west  profiles  the  com- 
bination of  4g  , T , , and  TiA  is  very  close  to  the  optimal 

r <j>  <p 
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Fig. 5.1  Local  covariance  function 

Standard  error  of  downward  continuation 
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solution.  However,  for  profile  spacings  smaller  than  about  0°.35 
the  combination  of  Ag  , Tr.  , and  T ^ gives  better  results. 

This  shows  that  the  high  resolution  of  the  T ^ ; -grad i ent s will 

only  be  effective  in  a limited  area  below  the  measuring  point. 

For  blocks  of  small  size  this  can  advantageously  be  used  because 
it  will  not  only  give  a better  accuracy  but  also  reduce  the  corre- 
lation between  neighbouring  points.  If  the  point  to  be  inter- 
polated is  directly  below  the  flight  path,  the  combination  of 
Ag  , T , and  T for  a wide  point  spacing  (>  0J.35)  and  of 

Ag  , T ^ , and  for  a dense  point  spacing  (>  0°.35)  will 

give  the  best  results.  For  north-south  profiles  T , T have 

r 'I1  $ $ 

to  be  replaced  by  T , T)(  and  vice  versa  to  get  the  corres- 
ponding results. 

Some  experiments  have  been  made  with  a different  arrange- 
ment of  profiles.  If  a constant  spacing  is  kept,  the  addition 
of  more  profiles  will  affect  the  results  only  for  spacings  below 
0°.3,  i.e.  if  the  ratio  between  profile  spacing  and  flying  alti- 
tude is  smaller  than  3:1.  Even  then,  the  effect  will  not  be  very 
marked,  so  that  the  results  presented  in  figures  5.1  and  5.2  are 
fairly  general.  The  use  of  a grid  of  east-west  and  north-south 
profiles  has  also  be  considered.  With  a fixed  number  of  flights 
the  spacing  of  the  profiles  for  the  grid  system  will  be  much 
wider  than  for  profiles  in  one  direction  and  the  interpolation 
accuracy  will  be  worse  in  this  case.  Therefore,  cross  profiles 
should  only  be  used  for  updating  purposes. 

For  further  reference  the  main  results  of  the  last  two 
sections  are  summarized: 

a. )  The  accuracy  of  an  interpolated  point  is  mainly  dependent  on 

its  distance  to  the  nearest  profile.  Therefore,  a dense 
profile  spacing  is  more  important  for  the  accuracy  of 
interpolation  than  cross  profiles. 

b. )  The  dense  sequence  of  points  along  the  profiles,  necessary 

when  integrating  first-order  gradients,  is  not  needed  for 
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c.)  For  profiles  with  arbitrary  azimuth  an  optimal  data  combi- 
nation should  include: 

sg  , T , , T , T , T for  profile  spacings  larger  than  0°.3 

r $ r a ' 

Ag  , T r , Trr  for  profile  spacings  smaller  than  0° . 3 . 

The  inclusion  of  more  second-order  gradients  will  only  slightly 
improve  the  accuracy  but  may  add  to  the  numerical  problems. 

A smaller  number  of  gradients  may  be  used  in  special  cases. 


6.  Influence  of  Measuring  Errors. 

It  has  been  mentioned  in  the  introduction  that  the  results 
of  this  accuracy  study  hinge  on  realistic  assumptions  about  the 
structure  of  the  gravity  field  and  the  accuracy  of  the  measure- 
ments. At  the  moment  no  airborne  gradiometer  is  in  operation 
and  therefore  realistic  estimates  are  difficult  to  obtain.  It 
is  expected  that  the  instruments  will  give  a standard  error  of 
+ 1 E for  any  component  integrated  over  an  interval  of  10  sec. 
Gravity  values  at  flight  level  can  be  obtained  by  integrating 
the  gradiometer  values  along  the  flight  path.  Moritz  (1975)  has 
studied  the  error  propagation  of  Ag  along  a profile  under 
certain  simplifying  assumptions.  The  standard  error  of  Ag  is 
dependent  on  the  length  of  the  profile  and  varies  between  + 1.4 
mgal  for  a distance  of  100  km  and  + 4.3  mgal  for  a distance 
of  1000  km  when  using  the  above  standard  error  for  the  gradio- 
meter  data.  Accurate  initial  values  and  a negligible  influence 
of  the  interaction  term  are  assumed  in  this  case.  Thus,  a value 
of  + 1 mgal  can  be  regarded  as  a fairly  optimistic  estimate  for 
the  standard  error  of  the  measurements.  Actual  errors  may  be 
considerably  larger. 

In  order  to  get  an  idea  about  the  sensitivity  of  the 
Ag-estimation  with  respect  to  errors  in  the  different  measure- 
ments, the  standard  errors  for  each  type  of  data  have  been  varied 
separately.  Again,  formula  (4.3)  can  be  used.  The  variances  o . 
in  equation  (4.4)  are  the  quantities  which  must  be  changed  in 
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this  case.  Fig.  6.1  shows  the  standard  error  of  estimation  m 

3 S 

as  a function  of  the  standard  error  m of  the  gravity  anomalies 

and  errorless  gradiometer  measurements.  Conversely,  errorless 

Ag-values  are  assumed  in  fig.  6.2  and  the  error  of  estimation  is 

a function  of  the  standard  error  m in  the  gradiometer  data. 

gg 

Results  are  shown  for  different  profile  spacings.  The  curves 

indicate  that  m is  more  sensitive  to  errors  in  the  second  - 

s 

order  gradients  than  to  errors  in  the  gravity  anomalies.  In 
fig.  6.1  the  curves  are  approximately  parallel  and  the  increase 
in  m is  less  than  half  of  that  in  the  measurements.  In  fig. 

S 

6.2  the  slope  of  the  curve  increases  sharply  with  growing 
profile  spacing.  At  a standard  error  = + 10  E the  curve 

levels  off  to  the  value  of  pure  Ag- i nter pol ati on  , i.e.  the 
gradiometer  data  are  not  accurate  enough  to  influence  the  error 
of  estimation.  Thus,  a high  accuracy  in  the  gradiometer  data  will 
control  even  a large  error  in  Ag  , while  large  errors  in  the 
second-order  gradients  usually  cannot  be  balanced  by  highly 
accurate  gravity  anomalies. 

The  situation  is  somewhat  different  for  small  profile 
spacings.  Here  the  increase  in  the  standard  error  of  estimation 
is  moderate  in  both  cases.  Thus,  a large  error  in  one  type  of 
data  is  balanced  by  a high  accuracy  of  the  other  one.  To  obtain 
optimal  results,  the  high  resolution  of  the  gradiometer  values 
for  small  distances  must  be  matched  by  a high  accuracy  of  the 
gravity  a noma  lies. 

The  difference  of  the  results  for  small  and  for  large 
profile  spacings  is  mainly  due  to  the  difference  in  the  corre- 
lation length  of  the  two  quantities  involved.  An  independent  check 
is  provided  by  results  published  in  Moritz  (1974).  Generally,  it 
can  be  concluded  that  for  profile  spacings  larger  than  0°.3  the 
accuracy  of  the  gradiometer  values  is  most  important  while  for 
profile  spacings  smaller  than  0°.3  a high  accuracy  of  both 
quantities  is  essential. 


I 
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7 . Mean  Values. 


So  far  we  have  only  considered  the  determi nation  of  point 
gravity  anomalies.  Usually  mean  values  for  blocks  of  specified 
size  are  required.  The  modifications  of  the  basic  formula  (4.3)  for 
this  case  are  discussed  in  this  section. 

Meissl  (1971)  has  shown  how  to  derive  a mean  anomaly  co- 
variance  function  from  a point  anomaly  covariance  function  by 
using  a narrow-sense  smoothing  operator.  Such  an  operator  is  de- 
fined by  the  properties  X ->-0  and  X = 1 , where  x.  are  the 

n o l 

eigenvalues  of  the  operator.  The  averaging  operator  over  a 
circular  cap  of  the  sphere  with  angular  radius  i//  is  an 
operator  of  this  kind.  Its  eigenvalues  are  given  by 
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or  by  the  recurrence  relation 
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6 = i - It-  tttt  , (cos*  ) - P . ( c o s ip  )}.  7.2 
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Replacing  a block  of  rectangular  form  by  a circular  cap  of  equal 
area  we  can  compute  a mean  anomaly  covariance  function  by  apply- 
ing the  smoothing  operator  to  the  point  covariance  function. 
Since  in  equation  (3.12)  we  have  the  spectral  representation  of 
the  point  anomaly  covariance  function  and  in  equation  (7.1)  the 
eigenvalues  of  the  smoothing  operator  we  can  directly  write 


K ( P , Q ) = A J e2k  sn  + 1P  ( co  sip ) 7.3 

u n ii  n 

11=  o 

where  K(P,Q)  is  the  mean  anomaly  covariance  function  derived 
from  K(P,Q).  Blocks  of  different  sizes  will  produce  different 


sets  of  a -values.  In  this  way,  we  can  check  once  more  the 

It 

model  covariance  function  against  empirical  values.  Taking  e.g. 
the  regional  covariance  function  with  all  coefficients  we  get 
a variance  of  C = 311  mgal*  for  5 x 5°blocks  and  a variance 

O 

of  C = 996  mqal"  for  l°x  1 blocks  as  compared  to  302  mgal2 

° 2 

and  996  mgal  from  empirical  data.  Although  such  a comparison 
is  not  very  sensitive  with  respect  to  the  local  behaviour  of  the 
covariance  function,  it  gives  a good  check  on  its  global  validity. 

For  all  practical  purposes  the  infinite  series  in  equation 
(7.3)  can  now  be  replaced  by  a finite  series.  For  a certain 
number  N the  terms  n > N will  not  contribute  significantly  to 
the  value  of  the  covariance  function.  Evidently,  the  smoothing 
property  will  be  stronger  for  large  blocks,  i.e.  the  number  11 
will  be  smaller  for  the  covariance  funct  on  of  5 x 5 blocks 
than  for  1°  x 1°  blocks.  Table  7.1  shows  the  number  of  necessary 
coefficients  for  different  block  sizes.  N has  been  determined 
such  that  the  terms  n > N will  not  change  the  integer  values  of 
the  covariance  function. 


Block  size 
N 


3000 

1200 

i 

400 

200 

Table  7.1  Number  of  necessary  coefficients  6 for 

n 

different  block  sizes. 

But  even  with  a limited  number  of  coefficients  the  actual 
summation  of  the  series  (7.3)  will  usually  cost  too  much  computer 
time.  It  would  be  advantageous  if  we  could  replace  the  infinite 
series  (7.3)  by  a closed  covariance  expression  as  it  can  be  done 
for  the  series  (4.3).  This  is  not  possible  with  the  B -values 
as  defined  by  formula  (7.1).  We  will  therefore  try  to  replace 
these  coefficients  by  quantities  which  make  such  a summation 
feasible  and  are  still  close  enough  to  the  exact  8^-values. 


Apparently  it  will  be  sufficient  to  have  a good  approxi- 
mation for  the  range  of  values  given  in  table  7.1  if  e -»  0 

II 

can  be  secured  for  n > N . Keeping  in  mind  that  the  g ar< 

II 

only  defined  for  integer  values  of  n , we  can  represent  them 
as  in  fig.  7.1  where  the  full  line  gives  the  6 -curve  for  the 
1°  x 1°  blocks. 


Fig. 7.1  Smoothing  coeff ici ents  6 and  approximations. 

It  is  quite  interesting  that  the  first  zero  will  roughly  coin- 
cide with  the  number  N given  in  table  7.1,  i.e.  that  the  oscilla 
ting  B -values  will  not  show  up  significantly  in  the  covariances 
This  is  true  for  all  block  sizes  smaller  than  2°  x 2°. For  5°  x 5° 
blocks  the  oscillating  part  contributes  about  1 %.  Thus,  for 
practical  purposes  it  will  be  sufficient  to  approximate  the 
g^ -curve  up  to  its  first  zero  by  some  other  smoothing  curve.  A 
very  good  approximation  of  the  6 -values  in  the  required  range 
can  be  achieved  by  the  function  sinx/x  where 
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and  N is  the  number  of  the  coefficient  closest  to  the  first  zero. 
In  that  case 


e 

II 


N si n(n • n/N) 
ii  • n 


7.4 


and  as  can  be  seen  from  the  dashed  line  in  fig.  7.1  there  is  an 

excellent  fit.  Unfortunately,  there  is  no  obvious  way  to  get 

these  6 -coefficients.  Therefore,  a 
11 

simple  a pprox i ma t i on  has  been  chosen. 


7.5 

the  smoothing  property.  Then 

s)"+1Pii(cos1|;)  7.6 

to  obtain  closed  expressions  in  the 
same  way  as  for  the  series  (4.3).  This  approach  gives  a simple 
interpretation  of  the  averaging  process.  Because  of  u < 1 
the  radius  of  the  Bjerhammar  sphere  is  diminished  to  obtain 
the  mean  anomaly  covariance  function.  This  will  result  in 

smaller  values  C and  G and  ,n  a larger  correlation  length 

o o 

5 . Similarly,  using  the  definition  of  s , equation  (7.6)  may 
be  explained  as  an  upward  continuation  of  the  original  covariance 
function  to  a certain  altitude  without  changing  the  radius  of  the 
Bjerhammar  sphere.  This  interpretation  has  been  used  by  Tscherning 
and  Rapp  (1974)  and  basically  their  method  to  obtain  a mean 
anomaly  covariance  function  is  equivalent  to  equation  (7.6). 
Because  of  this  dual  i nterpretat i on  it  can  be  expected  that  for 
a given  block  size  the  effects  of  downward  continuation  and  of 
mean  value  determination  will  cancel  at  a certain  altitude  h . 


closed  expressions  with 
somewhat  crude  but  very 
We  set 


2 n + 1 

6 = u 


where  u < 1 to  secure 


K(P,Q)  = A l k(u 


and  there  is  no  problem 


Table  7.2  shows  these  altitudes  for  different  block  sizes 


Block  size 


974153 


996234 


998884 

3.6 


999951 

0.2 


Table  7.2  Mean  values  and  levels  of  upward  continuation 


Take  as  an  example  the  estimation  of  mean  values  of  1°  x 1° 
blocks  at  sea  level.  At  a flight  level  of  12  km  the  effects 
of  downward  continuation  and  of  mean  value  determination  will 
cancel,  i.e.  they  can  be  replaced  by  planar  interpolation  at 
h = 12  km.  Usually  it  will  not  be  feasible  to  fix  flying  alti- 
tudes in  this  way  but  the  properties  of  a covariance  function  at 
a certain  altitude  may  be  appraised  by  comparing  it  to  the 
corresponding  mean  value  function.  Taking  as  before  h = 10  km, 
we  can  say  that  the  ag-values  at  this  level  will  be  almost  as 
smooth  as  1°  x 1°  mean  values  on  ground  level.  A very  dense 
spacing  of  points  will  therefore  increase  correlations  without 
essentially  improving  the  accuracy  of  the  results. 

Fig.  7.1  shows  two  8 -functions  of  type  (7.6)  as  full 
lines.  One  has  been  used  by  Tscherning  and  Rapp  (1974),  the  other 
one  corresponds  to  the  value  u given  in  table  7.2.  Both  curves 
are  only  very  crude  approximations  of  the  given  function  but  the 
values  of  the  corresponding  mean  anomaly  covariance  functions  are 
very  close  to  the  exact  ones  . From  the  difference  in  the  shape  of 
the  6 -curve  and  that  of  the  approximations  it  must  be  expected 
that  the  curvature  parameter  x will  be  changed  by  the  approxi- 
mation. Therefore,  the  resulting  covariance  function  should  not 
only  be  compared  for  Cq  but  at  least  up  to  a spherical  dis- 
tance of  1.5  • £ . This  has  been  done  for  the  8 -values  used  in 

n 

this  report  and  the  approximation  is  very  close  for  the  whole 
range.  Thus,  a more  refined  model  for  the  8 is  not  necessary 
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as  long  as  mean  gravity  anomalies  are  considered.  The  situation 
will  be  quite  different  for  mean  values  of  the  second-order 
gradients. 

In  order  to  apply  formula  (4.3)  we  still  need  the  cross- 
covariance function  between  point  and  mean  gravity  anomalies 
1<(P,Q).  In  analogy  to  formula  (7.3)  we  may  write 

CD 

K(P,Q)  = A l BiikuS,1+1Pn(cos*)  7.7 

n = o 


and  with  the  approximation  (7.5)  we  get 


oo 

I(P,Q)  = A l kfi(  /IT*  s)“+1  P (cos<|i)  . 7.8 

11=0 


We  then  have 


E = C - C C"lCT  7 . 9 

ss  ss  sx  sx 


for  the  estimation  of  mean  values.  Such  a formula  has  e.g.  been 
used  by  Lachapelle  (2975)  in  the  combination  of  gra v inet r i c and 
astrogeodetic  data.  It  should  be  noted,  however,  that  some  care 
must  be  taken  when  applying  it.  The  estimation  procedure  expressed 
by  formula  (7.9)  has  still  the  main  characteristics  of  an  inter- 
polation, i.e.  the  accuracy  is  dependent  on  the  distance  to  the 
nearest  profile.  Therefore,  the  mean  value  should  not  be  defined 
as  the  value  at  the  centre  of  the  block  because  otherwise  the 
accuracy  estimate  will  be  misleading  if  e.g.  a profile  passes 
through  the  centre.  Instead,  a regular  grid  of  points  covering 
the  area  of  the  block  can  be  used  to  determine  the  mean  standard 
error.  If  properly  chosen,  the  number  of  points  must  not  be  larger 
than  about  twice  the  number  of  profiles  passing  through  the  block. 
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Furthermore,  the  downward  continuation  problem  can  now  be  ex- 
plained in  a more  penetrating  way.  The  covariance  function  of  the 
measurements  at  flight  level  is  close  to  a mean  anomaly  co- 
variance  function  of  1°  x 1°  blocks.  If  we  want  to  estimate  mean 
values  of  smaller  block  sizes  or  even  point  values  from  these 
observations  we  must  be  aware  of  the  fact  that  we  are  trying  to 
reverse  the  averaging  process.  The  results  obtained  in  this  way 
will  be  very  smooth  and  strongly  correlated.  Thus,  we  cannot 
expect  to  recover  the  actual  point  values.  Furthermore,  it  is 
well  known  that  stability  problems  will  arise  in  such  unsmoothing 
procedures  and  that  measuring  errors  at  flight  level  may  be 
strongly  amplified,  see  Schwarz  (1973). 

Keeping  these  reservations  in  mind,  table  7.1  can  be  inter- 
preted. Mean  gravity  anomalies  are  given  for  different  block 
sizes  using  various  profile  spacings  and  different  assumptions 
on  the  measuring  errors.  The  best  combination  of  data  has  been 
determined  in  each  case  by  the  results  summarized  at  the  end  of 
section  5.  There  are  three  important  comclusions: 

a. )  The  accuracies  for  all  block  sizes  are  in  a range  which 

will  make  them  very  useful  for  geodetic  applications.  With 
a profile  spacing  of  1°  standard  errors  of  + 3 mgal  for 
5°  x 5°  blocks  and  of  _+  5 mgal  for  1°  x 1°  blocks  can  be 
achieved.  With  a profile  spacing  of  0°.3  standard  errors  of 
+ 3 mgal  for  15'  x 15'  and  5'  x 5‘  blocks  are  obtainable.  It 
should  be  noted,  however,  that  correlations  will  be  larger 
for  the  5'  x 5’  means  especially  if  the  standard  errors  of 
the  Ag-values  are  large. 

b. )  The  influence  of  measuring  errors  is  more  pronounced  for 

small  blocks  than  for  large  ones.  This  becomes  even  more 
apparent  if  always  the  same  data  combination  is  used  which 
has  not  been  done  in  table  7.1. 

c. )  The  difference  between  the  regional  and  the  local  covariance 

function  shows  clearly  for  large  block  sizes.  Standard  errors 
are  about  the  same  for  blocks  smaller  than  1°  x 1°.  For 
larger  blocks  the  results  obtained  by  the  regional  function 


Block  Profile  Regional  Covariance  Function  Local  Covariance  Function 

Size  Distance  Standard  Errors  in  mgal  Standard  Errors  in  nigal 


CO 

O 

CO 

CO 

C\J 

+ 1 

+ 1 

4-| 

lO  r-4  lO 

OO  C\J  «-H 


CO 

*— 4 

in 

in 

CO 

+ 1 

+ i 

+ 1 

oo 

1 — 4 

CO 

to 

CO 

+ 1 

+ i 

+ 1 

oo 

CO 

f— 4 

+ i 

+ i 

+ 1 

cn 

oo 

*3- 

co 

CO 

oo 

4-| 

+ 1 

+ 1 

to 

o 

cn 

OJ 

oo 

f—i 

+ 1 

+ 1 

4-| 

CM 

00 

LO 

t-H 

rH 

+ 1 

+ 1 

4*  1 

LO  O 
LO  CO 


oo  c\j  o 

00  to  LO 

+1  +1  4| 


O-v 

CO 

oo 

oo 

o 

i-H 

lO 

in 

co 

oo 

co 

CO 

in 

CO 

OJ 

CO 

oo 

+ 1 

+ 1 

+ 1 

+ 1 

4-1 

+ 1 

4-1 

4-1 

4-1 

oo 

CO 

CO 

+i 

+i 

+ 1 

o°. 

in 

in 

o . 

oo 

o 
0 • 
r—1 

O lO 

o • o • 

•“»  O 


r-, 

r-^ 

OO 

ro 

T—i 

r— 4 

+ 1 

4-1 

4-1 

in 

co 

OO 

O' 

o- 

o • 

o 

o 

o 

in 

f—* 

43 


are  certainly  more  reliable.  If  we  assume  isostatic  anomalies 
for  the  sea  surface  their  covariances  are  probably  well  re- 
presented by  the  local  covariance  function. 

8.  Combination  with  Satellite  Altimetry. 

So  far,  we  have  estimated  ig-values  from  first  and  second 
order  gravitational  gradients  only.  When  using  measurements  from 
satellite  altimetry  we  are  not  just  adding  a new  set  of  data  but 
we  are  introducing  a very  different  type  of  observations.  They 
belong  to  a class  of  functions  which  is  smoother  than  that  of 
the  estimated  values.  This  can  simply  be  seen  from  Stokes’ 
integral  where  geoidal  undulations  are  determined  by  integrating 
gravity  anomalies.  Reversing  the  process  of  integration  will 
lead  to  a number  of  intricate  mathematical  problems.  It  is  this 
problem  we  are  dealing  with  and  the  difficulties  connected  to 
the  so-called  inversion  of  Stokes  integral  are  hidden  in  our 
approach,  too. 

The  smoothness  of  the  altimeter  data  will  produce  a sort 
of  basic  standard  error  in  the  A g -est i ma t i on . By  this  we  mean  that 
even  with  an  ideal  data  distribution  and  errorless  data  we  will 
not  be  able  to  push  the  standard  error  of  estimation  below  a 
certain  limit.  Furthermore, numeri cal  instabilities  will  occur 
if  the  point  density  exceeds  a certain  measure.  From  practical 
considerations  it  is  important  to  find  the  optimal  point  dis- 
tance where  the  basic  standard  error  can  be  obtained  with  a 
minimum  number  of  points  and  profiles.  Also,  the  critical  dis- 
tance where  instabilities  are  likely  to  occur  should  be  known 
in  order  to  provide  safeguards  in  the  program. 

The  same  configuration  as  in  section  4 and  5 has  been  used 
to  estimate  gravity  anomalies  from  geoidal  undulations  and  to 
determine  the  two  c hara cter i st i c distances.  The  standard  error 
of  the  altimeter  data  has  been  varied  between  0 m and  + 5 m. 
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Only  the  regional  covariance  function  has  been  used  in  this 
case.  Results  are  shown  in  fig.  8.1  where  the  error  of  esti- 
mation is  given  as  a function  of  the  profile  distance  for 
different  accuracies  of  the  altimeter  data.  The  numbers  written 
to  the  curves  are  standard  errors  in  meter  units.  Let  us  first 
consider  pure  interpolation,  i.e.  the  curve  numbered  0.  The 
standard  error  levels  off  at  a profile  distance  of  about  0?2 
and  becomes  erratic  at  about  0?08.  Profile  distances  smaller 
than  this  critical  value  will  give  completely  irregular  results 
including  imaginary  standard  errors.  The  critical  value  is 
dependent  on  the  internal  accuracy  of  the  computer  and  there  is 
also  some  dependence  on  the  point  configuration.  Thus,  for  double 
precision  arithmetic  a critical  distance  of  0? 1 is  necessary  to 
secure  a stable  solution.  On  the  other  hand,  the  optimal  distance 
is  at  about  0.2  and  measurements  with  a denser  profile  spacing 
will  not  improve  the  solution.  For  mean  values  the  critical 
distance  stays  the  same,  while  the  optimal  distance  increases 
with  growing  blocks  size.  Thus,  we  have  optimal  distances  of 
about  0:4  and  0^3  for  l°x  1J  and  5°x  5°  blocks.  Using  these 
values  we  are  able  to  determine  the  basic  standard  error  which 
is  + 19.5  mgal  for  point  values,  + 6.5  mgal  for  l°x  1°  blocks, 
and  + 3.0  mgal  for  5°x  5°  blocks.  For  blocks  smaller  than  l°x  lc, 
correlations  become  so  large  that  it  is  misleading  to  give  stan- 
dard errors  as  estimates  of  accuracy. 

If  we  turn  from  pure  interpolation  to  interpolation  inclu- 
ding measuring  errors  we  can  observe  some  interesting  features 
from  fig.  8.1.  Growing  standard  errors  of  the  altimeter  data  will 
result  in  a larger  optimal  distance  and  in  a flattening  of  the 
error  curves.  At  m = + lm  we  have  an  almost  linear  slope  for 
the  error  of  estimation.  -Judging  from  the  set  of  curves  the 
optimal  distance  is  at  about  1°  in  that  case  and  the  apparent 
increase  in  accuracy  with  smaller  profile  distances  must  be  re- 
garded as  fictitious.  Furthermore,  we  have  stable  solutions  even 
for  very  small  profile  distances.  Thus,  the  incorporation 


ttk 


Fig. 8.1  Estimation  of  Ag-values  from  altimeter  data  only 


4 


of  measuring  errors  in  the  equation  system  will  stabilize  the 

solution  by  smoothing  it.  Only  for  standard  errors  of  about 

m = + 0.1  m we  can  expect  a solution  which  will  recover  some 
n - 

detail,  an  accuracy  of  m = + 1 m will  already  produce  a strongly 

N — 

smoothed  solution.  Similar  results  are  obtained  for  the  1° x 1° 
mean  anomaly  function,  except  that  the  curves  are  flatter  on  the 
whole.  For  5°x  5°  means  the  curves  are  already  very  smooth  and 
the  difference  between  the  +_  0.1  m and  the  + 1 m curve  is  not 
very  pronounced.  This  could  be  expected  because  in  case  of  large 
blocks  the  basic  standard  error  will  dominate  the  shape  of  the 
error  curve.  The  standard  errors  of  estimation  for  different 
blocks  sizes  and  different  accuracies  of  the  input  data  are  given 
in  table  8.1: 


Block  Standard  error  of  estimation  in  mgal  for  m = 


+ 0 . 1 m + 1 . 0 m 


+ 5.0m 


5°  x 5° 
l°x  1° 


+ 10.0  m 


+ 7.9 


+ 24.0 


Table  8.1:  Standard  errors  of  Ag-estimation  from  altimeter 
data  only. 


The  results  are  somewhat  different  from  those  published  by  Rapp 
(1974).  This  is  probably  due  to  differences  in  the  covariance 
functions. 

So  far,  we  have  considered  altimeter  data  as  the  only 
available  data  set  for  the  estimation  of  Ag  . Let  us  now  turn 
to  combination  solutions.  The  results  on  the  critical  and  optimal 
distance  will  remain  valid  for  this  case,  too.  But  they  have  to 
be  supplemented  using  a somewhat  different  viewpoint.  We  are  now 
asking  for  that  profile  distance  where  altimeter  data  cease  to 
contribute  significantly  to  the  solution.  Again,  this  question 
is  of  practical  importance  because  it  will  help  to  reduce  the 
number  of  observations  and  to  avoid  unnecessary  correlations. 


' 
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Fig.  8.2  shows  the  error  of  estimation  as  a function  of 
the  profile  distance  for  combination  solutions  with  and  without 
altimeter  data  and  for  blocks  of  different  sizes.  The  standard 
errors  of  the  altimeter  measurements  have  been  varied  while 
standard  errors  of  the  other  quantities  have  been  held  fixed. 

Point  estimation  has  been  included  as  a limiting  case  and  because 
it  permits  a direct  comparison  with  previous  results.  Since  the 
simple  configuration  as  in  fig.  8.1  has  been  used  this  is  not  true 
for  the  mean  values.  Their  standard  errors  will  need  a scaling. 

The  standard  errors  of  the  altimeter  data  have  again  been  written 
to  the  curves.  There  are  three  obvious  results.  First,  the  accu- 
racy of  the  5x  5°  mean  values  will  always  be  improved  when  using 
satellite  altimetry.  Second,  for  blocks  of  l°x  1°  and  smaller, 
satellite  altimetry  will  not  contribute  significantly  if  the  pro- 
file spacing  is  smaller  than  1°  . It  is  also  at  this  distance 
that  the  standard  error  of  the  5°x  5°  blocks  levels  off  indi- 
cating that  a denser  spacing  will  not  improve  the  accuracy.  Third, 
an  improvement  in  the  altimeter  accuracy  from  mN  = + 1.0  m to 
m^  = 0.1  m will  improve  the  accuracy  of  the  combination  solu- 

tion by  only  about  5 % . These  results  have  been  checked  and 
confirmed  by  computing  a number  of  cases  in  table  7.1  from  a 
combination  of  gravity,  gradiometer,  and  altimeter  data.  Finally, 
it  should  be  remembered  that  the  standard  errors  of  Ag  and  of 
the  second-order  gradients  have  been  chosen  rather  large.  If  they 
are  smaller,  the  contribution  of  satellite  altimetry  to  blocks  of 
l°x  1°  and  smaller  sizes  will  be  even  less  significant. 


9.  Numerical  Problems. 


In  this  section  we  will  first  treat  the  question  of 
numerical  stability  in  a more  comprehensive  way  and  then  discuss 
some  simplifications  in  the  inversion  of  the  covariance 
matrix  . Again,  questions  of  practical  solvability  will  be 

more  important  than  mathematical  rigour. 
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Instabilities  nay  occur  in  two  cases.  First  in  the  downward 
continuation  from  flight  level  to  sea  level  and  second  in  the 
estimation  of  Ag  from  altimeter  data  only.  In  the  first  case 
the  smoothness  of  the  data  function  is  due  to  attenuation  and 
dependent  on  the  size  of  s in  equation  (3.12).  In  the  second 
case  the  data  function  has  a lower  order  of  smoothness  than 
the  runction  to  be  estimated.  This  is  apparent  from  the  coef- 
ficients k;  in  the  covariance  models.  From  a practical  point 
of  view  both  problems  can  be  treated  in  a similar  way  because 
the  quantity 


I 


!i 


t 
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9.1 


represents  a damping  of  the  original  coefficients  in  both  models 
and  therefore  high  frequencies  cannot  properly  be  recovered. 

If  we  have  an  unstable  linear  system 

x = As  9.2 

we  can  regularize  the  solution  by  augmenting  the  diagonal  of  the 
system  of  normal  equations  by  small  amounts,  i.e.  by  adding  a 
matrix  al  . We  obtain 


s = (AtA  + a I ) - 1 ATx 


9.3 


as  a solution  and  the  usual  minimum  condition  of  least-squares 
adjustment  is  replaced  by 

||  A s - x || 2 + a • ||  s || 2 -*■  m i n i mum  9.4 

where  ||  • || 2 denotes  the  square  norm.  This  method  of  regulari- 
zation is  due  to  Tihonov  (1965)  and  it  is  valid  ev  n in  cases 
when  the  operator  A is  nonlinear.  As  can  be  seen  from  the 
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condition  (9.4)  strong  variations  of  the  solution  are  prevented 
by  imposing  a restriction  on  the  norm  of  the  solution.  In  practice, 
the  determination  of  a is  a difficult  problem,  see  Ivanow  (1966). 
If  it  is  chosen  too  small,  we  will  have  residual  instabilities. 

If  it  is  made  too  large,  the  solution  will  be  too  smooth  and  the 
optimum  will  not  be  reached.  Usually  one  has  to  rely  on  trial  and 
error.  An  interesting  numerical  example  is  given  in  Tikhonov  and 
Glasko  (1964)  where  some  of  the  computational  difficulties  are 
di scussed  . 

It  has  been  pointed  out  in  Sc hwarz  (1973)  that  in  case  of 
a linear  system  there  is  a close  connection  between  Tihonov's 
method  and  least-squares  collocation.  A similar  argument  can  be 
used  in  our  case.  The  variances  c.  in  equation  (4.4)  can  be 
interpreted  as  regularizing  factors.  However,  they  are  given  in 
advance  by  the  accuracy  of  the  data  and  have  not  been  determined 
by  a best  approximation  of  some  kind.  It  is  therefore  still 
necessary  to  find  an  optimal  coefficient  a and  to  see  by  com- 
parison how  much  we  loose  by  being  bound  to  the  given  c . In 
other  words,  because  we  must  use  an  actual  variance  our  solution 
will  usually  be  too  smooth.  Only  if  our  data  are  given  with  a 
very  high  accuracy  we  can  hope  to  approach  the  optimal  a . To 
give  more  practical  meaning  to  this  concept  let  us  return  to 
fig.  8.1.  The  coefficients  a can  be  obtained  by  squaring  the 
standard  errors  and  it  is  quite  obvious  that  an  optimal  solution 
must  have  a<0.01  . The  loss  of  information  with  growing  a is 
considerable  and  apparently  the  difference  between  a=0.01  and 
a= 1 . 0 is  not  only  of  a quantitative  but  also  of  a qualitative 
manner.  More  detailed  results  can  be  obtained  from  simulation 
studies  where  the  exact  functions  are  known  and  an  independent 
control  is  possible. 

A different  approach  to  the  problem  is  provided  by  looking 
at  its  spectral  representation.  Using  formulas  (3.12)  and  (9.1) 
we  can  write 


K(P.Q)  = A l t^PJcos*) 
n*o 


9.5 
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where  the  coefficients  t are  a smoothed  version 

n 

k;  . The  smoothing  is  dependent  on  the  degree  n . 
by  a factor  v“  + 1 in  case  of  attenuation  and  by  a 
R /(n-1)  in  case  of  altimeter  data.  The  quantity 
by 


(R  + h)2 


where  h is  theflying  altitude.  In  our  estimation  procedure  we 
are  basically  trying  to  recover  the  original  coefficients  . 

Errors  in  the  data  will  falsify  the  coefficients  t by  small 

II 

amounts.  Let  us  assume  for  a moment  that  we  have  a small  error 
of  equal  size  in  each  coeficient  t . When  we  try  to  recover 
the  quantities  k , the  effect  of  this  error  will  be  very 

n 

different  depending  on  the  degree  of  the  coefficient.  For  low 
degrees  we  will  have  an  amplifying  factor  close  to  one  and  no 
difficulties  will  arise.  For  high  degrees  the  amplifying  factor 
may  be  arbitrarily  large.  Thus,  our  solution  will  become  completely 
meaningless  if  we  try  to  include  all  high  frequencies.  To  over- 
come this  difficulty  some  sort  of  filtering  must  take  place.  This 
filtering  procedure  is  the  analogue  to  the  determination  of  an 
optimal  a . In  one  case  we  are  considering  the  frequency  domain, 
in  the  other  case  the  time  domain.  The  filtering  of  the  high 
frequencies  will  as  before  produce  a smoothing  of  the  solution. 

The  extent  of  the  smoothing  will  depend  on  the  characteristics  of 
the  smoothing  factors  as  well  as  on  the  accuracy  level  of  the 
data.  In  our  problem  the  smoothing  of  the  high  frequencies  is  much 
stronger  in  case  of  satellite  altimetry.  The  accuracy  level  of  the 
data  can  be  related  to  the  effects  of  the  filtering  process  in  the 
frequency  doma in.  Eigenvalue  techniques  as  in  Schwarz  (1975)  could 
be  used  for  this  purpose.  They  will  allow  interesting  comparisons 
with  the  results  presented  in  fig.  8.1.  Thus,  it  is  important  to 
take  into  account  both  aspects  of  the  stability  problem  in  order 
to  give  a meaningful  interpretation  to  the  numerical  results. 


of  the  original 
It  is  generated 
factor 

v is  def i ned 


51 


Let  us  finally  consider  the  inversion  of  the  matrix  C 
in  equation  (4.1).  Since  the  inversion  is  the  time-consuming 
part  in  the  algorithm,  it  will  strongly  influence  the  efficiency 
of  the  program.  Some  general  rules  have  already  been  given  in 
section  4 how  to  reduce  the  number  of  measurements  and  by  this 
the  size  of  the  matrix.  We  shall  therefore  concentrate  on 
simplifications  arising  from  the  structure  of  the  matrix.  Table 
9.1  shows  the  covariance  matrix  for  the  optimal  data  combination 
as  given  in  section  5 including  altimeter  data.  Only  one  measured 
point  is  considered . 


The  variances  and  covariances  of  the  upper  half  of  the  symmetrical 
matrix  are  shown.  The  notations  are  as  follows: 

0 ...  covariance  is  zero  for  all  ijj 
(0)...  covariance  is  zero  for  = 0 
+ ...  variance  and  covariance  have  a value. 


The  crosscovariances  implying  T , ^ form  a diagonal  submatrix  which 


can  directly  be  inverted.  In  this  way  more  than  20  % of  the  com- 
puting time  for  the  inversion  may  be  saved.  Furthermore,  there 


are  some  zero-bands  connected  with  T and  zeros  from  other 

f r 


crosscovariances  may  be  matched  with  certain  profile  spacings. 
This  will  not  influence  the  speed  of  inversion  but  will  make  the 
procedure  more  stable  because  a number  of  off-diagonal  terms  will 
become  zero.  It  may  also  be  possible  to  introduce  in  advance  a 


- — - - — 


simple  form  of  pivoting  by  considering  the  size  and  the  shape  of 
the  different  covariance  functions  and  by  ordering  them  accordingly. 
Thus,  a careful  assessment  of  the  matrix  structure  may  help  to 
speed  up  and  stabilize  the  inversion. 


10 . Concl usions 


The  preceding  i nvestigations  have  shown  that  a combination 
of  first  and  second  order  gravitational  gradients  measured  at 
flying  altitudes  will  supply  mean  gravity  anomalies  at  ground 
level  accurate  enough  for  geodetic  applications.  The  results 
fall  roughly  into  two  groups.  Those  which  are  basically  inde- 
pendent of  the  chosen  covariance  function  and  those  which  are 
not.  Results  of  the  first  group  may  be  summarized  in  six  points: 

1.  The  accuracy  of  an  interpolated  gravity  anomaly  is  mainly 
dependent  on  its  distance  to  the  nearest  profile.  Therefore, 
a dense  profile  spacing  is  more  important  for  interpolation 
than  cross  profiles.  Cross  profiles  should  be  used  for  up- 
dating only. 

2.  A dense  sequence  of  points  along  the  profile  is  only  necessary 
for  integrating  first-order  gradients.  For  interpolation  of 
Ag-values  a density  between  1/4  and  1/2  of  the  profile  spacing 
will  be  sufficient. 

3.  For  profiles  with  an  arbitrary  azimuth  an  optimal  data  combi- 
nation should  include: 


Ag’  Tr4>  ’ TrX  ’ TH’ 

for 

1 a r ge 

prof i 1 e 

spacings  (> 

0°3) 

Ag  * Tr4>  ’ T r X * Trr 

for 

sma  1 1 

prof i 1 e 

spacings  (< 

0?3) 

The  inclusion  of  more  second-order  gradients  will  only  slightly 
improve  the  accuracy.  A smaller  number  of  gradients  may  be 
used  in  special  cases. 


4.  The  influence  of  measuring  errors  is  also  different  for  large 
and  for  small  profile  spacings.  In  the  first  case  the  accuracy 
of  the  gradiometer  values  is  most  important  while  in  the  second 
case  a high  accuracy  of  gravity  anomalies  and  second-order 
gradients  is  essential. 

5.  When  estimating  mean  values  the  covariance  function  has  to  be 
changed  in  a similar  way  as  for  upward  continuation.  Thus,  for 
a given  block  size  the  effects  of  downward  continuation  and  of 
mean  value  determination  will  cancel  at  a certain  altitude. 

6.  The  combination  of  first  and  second  order  gradients  with 
satellite  altimetry  will  essentially  contribute  to  the  estima- 
tion of  5°x  5°  means.  With  present  day  accuracy  only  marginal 
improvements  can  be  expected  for  l°x  1°  blocks.  For  blocks  of 
smaller  size  altimetry  should  not  be  used  because  correlations 
introduced  by  these  data  may  impair  the  accuracy. 

The  following  results  are  dependent  on  assumptions  made 
about  the  covariance  functions.  These  functions  have  been  varied 
in  a rather  wide  range  so  that  it  should  be  possible  to  estimate 
the  effect  of  a change  in  one  of  the  essential  parameters.  When 
actual  gradiometer  measurements  become  available,  the  basic 
assumptions  can  be  checked  and  appropriate  changes  be  made  if 
necessary.  Linder  these  premises  we  can  say: 

7.  Mean  gravity  anomalies  can  be  estimated  from  first  and  second 
order  gradients  at  a flight  level  of  10  km  with  the  following 
accuracies.  With  a profile  spacing  of  1°  standard  errors  of 

+ 3 mgal  for  5°x  5°  blocks  and  of  +_  5 mgal  for  l°x  1°  blocks 
can  be  achieved.  With  a profile  spacing  of  0?3  standard  errors 
of  + 3 mgal  are  obtainable  for  blocks  of  1 5 ' x 15'  and  5'x  5'. 

It  should  be  noted,  however,  that  correlations  may  be  large 
for  the  5'x  5'  blocks  depending  on  the  accuracy  of  the  Ag- 
values.  The  influence  of  measuring  errors  is  more  pronounced 
for  small  blocks  than  for  large  ones.  More  details  can  be 
found  in  table  7.1  of  this  report. 

The  method  used  in  this  accuracy  study  is  least-squares 
collocation.  Its  great  flexibility  has  been  demonstrated  in 
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handling  problems  as  different  as  the  combination  of  heterogeneous 
data,  interpolation,  downward  continuation  and  mean  value  esti- 
mation. Its  numerical  advantages  have  been  shown  in  discussing 
the  stability  problem  which  is  essential  in  the  treatment  of  the 
above  applications. 
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